
function [Wf,Wi,Wn,iterVF]=solve_value_function_A1(Par1,Ff_2,wf_2,Fi_2,wi_2,n)

r=0.036;

% Clenshaw-Curtis quadrature
% n+1 nodes: x in [-1,1] by decreasing order
[x,wcc] = ccquad(n-1); x=flipud(x);

transition_rates=Par1(1:6)
theta=Par1(11)
a=Par1(12)
gamma=0; % gamma = 0 at benchmark
b2 = Par1(13)
df_2=transition_rates(1);
di_2=transition_rates(2);
lnf_2=transition_rates(3);
lni_2=transition_rates(4);
lfi_2=transition_rates(5);
lif_2=transition_rates(6);


%preallocation to speed up

Wi = zeros(n,1); Rw_f_i = zeros(1,n); Rw_i_f = zeros(1,n); Wf_Rw_f_i = zeros(1,n); Wi_Rw_i_f = zeros(1,n);
Wf = zeros(n,1); Fi_2_Rw_i_f = zeros(1,n); Ff_2_Rw_f_i = zeros (1,n); derWf = zeros(1,n);derWi = zeros(1,n);
intmaxWf_Wi = zeros(1,n); intmaxWi_Wf = zeros(1,n); eWi = zeros(1,n); eWf = zeros(1,n);

%% Transforming the probabilities g and f in densities
% If wages are linearly spaced
wht_f_2=(wf_2(n)-wf_2(1))/(n-1); 
wht_i_2=(wi_2(n)-wi_2(1))/(n-1);


    % Initializing
    
    % value function (min and max)
    Wn=(1-exp(-theta*(min(wi_2))))/theta/r;
    Wf_min= (1-exp(-theta*(min(wf_2))))/theta/r; Wf_max= (1-exp(-theta*(max(wf_2))))/theta/r;
    Wi_min=(1-exp(-theta*(min(wi_2))))/theta/r; Wi_max=(1-exp(-theta*(max(wi_2))))/theta/r ;
    
     
    Wf=(Wf_min+Wf_max)/2+x*(Wf_max-Wf_min)/2; % Interpolation with CC quadrature (defined above)
    Wi=(Wi_min+Wi_max)/2+x*(Wi_max-Wi_min)/2;

    
    critVF=1;
    iterVF=1;
    omegaW=0.1;
    while mean(critVF)>0.005 && iterVF<=400
        
        % Reservation wages
        for i=1:n
            % to construct Wfn
          
            dWi_Wf=abs(Wi-Wf(i))./Wf(i);
            Rw_i_f(i)=wi_2(min(find(dWi_Wf<=min(dWi_Wf))));
            Wi_Rw_i_f(i)=Wi(min(find(dWi_Wf<=min(dWi_Wf))));
            Fi_2_Rw_i_f(i)=Fi_2(min(find(dWi_Wf<=min(dWi_Wf))));
        
            
            derWf(i)=1./(r+df_2 +lfi_2*Fi_2_Rw_i_f(i));
            
            % to construct Win
            dWf_Wi=abs(Wf(i,:)-Wi(i))./Wi(i);
            Rw_f_i(i)=wf_2(find(dWf_Wi<=min(dWf_Wi), 1 ));
            Wf_Rw_f_i(i)=Wf(find(dWf_Wi<=min(dWf_Wi), 1 ));
            Ff_2_Rw_f_i(i)=Ff_2(find(dWf_Wi<=min(dWf_Wi), 1 ));
            
            
            derWi(i)=1./(r+di_2+lif_2*Ff_2_Rw_f_i(i));
        end
       
                % to construct Wnn
                dWf_Wn=abs(Wf-Wn)./Wn;
                Rw_f_n=wf_2(find(dWf_Wn<=min(dWf_Wn), 1 ));
                Ff_2_Rw_f_n=Ff_2(min(find(dWf_Wn<=min(dWf_Wn))));
                
                dWi_Wn=abs(Wi-Wn)./Wn;
                Rw_i_n=wi_2(find(dWi_Wn<=min(dWi_Wn), 1 ));
                Fi_2_Rw_i_n =Fi_2(min(find(dWi_Wn<=min(dWi_Wn))));
        
        % solve integral by parts
        
        % Wf
        for i=1:n
            intmaxWi_Wf(i)=sum((wi_2>Rw_i_f(i)).*Fi_2.*derWi'.*wht_i_2); 
        end
        % Wi
        for i=1:n
            intmaxWf_Wi(i)=sum((wf_2>Rw_f_i(i)).*Ff_2.*derWf'.*wht_f_2);   
        end
        
        % Wnn
        
        intmaxWf_Wn=sum((wf_2>Rw_f_n).*Ff_2.*derWf'.*wht_f_2);
        intmaxWi_Wn=sum((wi_2>Rw_i_n).*Fi_2.*derWi'.*wht_i_2);
        
        %%%%% Value functions
        % solve in terms of F, reservation wages and lambdas
        
        % Wfn(w1) and Win(w1)
        for i=1:n
            
            eWf(i)=(1/(r+df_2))*((1-exp(-theta*(wf_2(i))))/theta+a+df_2*Wn ...
                +lfi_2*intmaxWi_Wf(i))-Wf(i);
            
            eWi(i)=(1/(r+di_2))*((1-exp(-theta*(wi_2(i))))/theta+gamma+di_2*Wn ...
                +lif_2*intmaxWf_Wi(i))- Wi(i);

            
        end

        % Wnn
        
        eWn=(1/r)*((1-exp(-theta*b2)))/theta+gamma ...
            +lnf_2*intmaxWf_Wn ...
            +lni_2*intmaxWi_Wn;

        
        % Update value functions
        
        Wf=Wf+omegaW*eWf';
        Wi=Wi+omegaW*eWi';
        Wn=Wi(1);eWn=0;
        
        critVF=[abs(eWf'./Wf);abs(eWi'./Wi);abs(eWn/Wn)];
        
        iterVF=iterVF+1;
     
    end
    iterVF
    Rw_2_f_ind = Rw_f_n; 
    Rw_2_i_ind = Rw_i_n;  
    save Rw_2_f_ind Rw_2_f_ind
    save Rw_2_i_ind Rw_2_i_ind
end
 